home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / SeqPups / appsrc / cap.src / CAP2.c next >
Text File  |  1996-07-05  |  63KB  |  2,421 lines
  1.  
  2. /* dgg -- for use as Mac child-app */
  3. #define MACINTOSH 1
  4.  
  5. /* CONTIG ASSEMBLY PROGRAM (CAP)
  6.  
  7.    copyright (c) 1991    Xiaoqiu Huang
  8.    The distribution of the program is granted provided no charge
  9.    is made and the copyright notice is included.
  10.  
  11.    Proper attribution of the author as the source of the software
  12.    would be appreciated:
  13.    "A Contig Assembly Program Based on Sensitive Detection of
  14.    Fragment Overlaps" (submitted to Genomics, 1991)
  15.     Xiaoqiu Huang
  16.     Department of Computer Science
  17.     Michigan Technological University
  18.     Houghton, MI 49931
  19.         E-mail: huang@cs.mtu.edu
  20.  
  21.    The CAP program uses a dynamic programming algorithm to compute
  22.    the maximal-scoring overlapping alignment between two fragments.
  23.    Fragments in random orientations are assembled into contigs by a
  24.    greedy approach in order of the overlap scores. CAP is efficient
  25.    in computer memory: a large number of arbitrarily long fragments
  26.    can be assembled. The time requirement is acceptable; for example,
  27.    CAP took 4 hours to assemble 1015 fragments of a total of 252 kb
  28.    nucleotides on a Sun SPARCstation SLC. The program is written in C
  29.    and runs on Sun workstations.
  30.  
  31.    Below is a description of the parameters in the #define section of CAP.
  32.    Two specially chosen sets of substitution scores and indel penalties
  33.    are used by the dynamic programming algorithm: heavy set for regions
  34.    of low sequencing error rates and light set for fragment ends of high
  35.    sequencing error rates. (Use integers only.)
  36.     Heavy set:             Light set:
  37.  
  38.     MATCH     =  2             MATCH     =  2
  39.     MISMAT    = -6             LTMISM    = -3
  40.     EXTEND    =  4             LTEXTEN   =  2
  41.  
  42.     In the initial assembly, any overlap must be of length at least OVERLEN,
  43.     and any overlap/containment must be of identity percentage at least
  44.     PERCENT. After the initial assembly, the program attempts to join
  45.     contigs together using weak overlaps. Two contigs are merged if the
  46.     score of the overlapping alignment is at least CUTOFF. The value for
  47.     CUTOFF is chosen according to the value for MATCH.
  48.  
  49.     DELTA is a parameter in necessary conditions for overlap/containment.
  50.     Those conditions are used to quickly reject pairs of fragments that
  51.     could not possibly have an overlap/containment relationship.
  52.     The dynamic programming algorithm is only applied to pairs of fragments
  53.     that pass the screening. A large value for DELTA means stringent
  54.     conditions, where the value for DELTA is a real number at least 8.0.
  55.  
  56.     POS5 and POS3 are fragment positions such that the 5' end between base 1
  57.     and base POS5, and the 3' end after base POS3 are of high sequencing
  58.     error rates, say more than 5%. For mismatches and indels occurring in
  59.     the two ends, light penalties are used.
  60.  
  61.     A file of input fragments looks like:
  62. >G019uabh
  63. ATACATCATAACACTACTTCCTACCCATAAGCTCCTTTTAACTTGTTAAA
  64. GTCTTGCTTGAATTAAAGACTTGTTTAAACACAAAAATTTAGAGTTTTAC
  65. TCAACAAAAGTGATTGATTGATTGATTGATTGATTGATGGTTTACAGTAG
  66. GACTTCATTCTAGTCATTATAGCTGCTGGCAGTATAACTGGCCAGCCTTT
  67. AATACATTGCTGCTTAGAGTCAAAGCATGTACTTAGAGTTGGTATGATTT
  68. ATCTTTTTGGTCTTCTATAGCCTCCTTCCCCATCCCCATCAGTCTTAATC
  69. AGTCTTGTTACGTTATGACTAATCTTTGGGGATTGTGCAGAATGTTATTT
  70. TAGATAAGCAAAACGAGCAAAATGGGGAGTTACTTATATTTCTTTAAAGC
  71. >G028uaah
  72. CATAAGCTCCTTTTAACTTGTTAAAGTCTTGCTTGAATTAAAGACTTGTT
  73. TAAACACAAAATTTAGACTTTTACTCAACAAAAGTGATTGATTGATTGAT
  74. TGATTGATTGATGGTTTACAGTAGGACTTCATTCTAGTCATTATAGCTGC
  75. TGGCAGTATAACTGGCCAGCCTTTAATACATTGCTGCTTAGAGTCAAAGC
  76. ATGTACTTAGAGTTGGTATGATTTATCTTTTTGGTCTTCTATAGCCTCCT
  77. TCCCCATCCCATCAGTCT
  78. >G022uabh
  79. TATTTTAGAGACCCAAGTTTTTGACCTTTTCCATGTTTACATCAATCCTG
  80. TAGGTGATTGGGCAGCCATTTAAGTATTATTATAGACATTTTCACTATCC
  81. CATTAAAACCCTTTATGCCCATACATCATAACACTACTTCCTACCCATAA
  82. GCTCCTTTTAACTTGTTAAAGTCTTGCTTGAATTAAAGACTTGTTTAAAC
  83. ACAAAATTTAGACTTTTACTCAACAAAAGTGATTGATTGATTGATTGATT
  84. GATTGAT
  85. >G023uabh
  86. AATAAATACCAAAAAAATAGTATATCTACATAGAATTTCACATAAAATAA
  87. ACTGTTTTCTATGTGAAAATTAACCTAAAAATATGCTTTGCTTATGTTTA
  88. AGATGTCATGCTTTTTATCAGTTGAGGAGTTCAGCTTAATAATCCTCTAC
  89. GATCTTAAACAAATAGGAAAAAAACTAAAAGTAGAAAATGGAAATAAAAT
  90. GTCAAAGCATTTCTACCACTCAGAATTGATCTTATAACATGAAATGCTTT
  91. TTAAAAGAAAATATTAAAGTTAAACTCCCCTATTTTGCTCGTTTTTGCTT
  92. ATCTAAAATACATTCTGCACAATCCCCAAAGATTGATCATACGTTAC
  93. >G006uaah
  94. ACATAAAATAAACTGTTTTCTATGTGAAAATTAACCTANNATATGCTTTG
  95. CTTATGTTTAAGATGTCATGCTTTTTATCAGTTGAGGAGTTCAGCTTAAT
  96. AATCCTCTAAGATCTTAAACAAATAGGAAAAAAACTAAAAGTAGAAAATG
  97. GAAATAAAATGTCAAAGCATTTCTACCACTCAGAATTGATCTTATAACAT
  98. GAAATGCTTTTTAAAAGAAAATATTAAAGTTAAACTCCCC
  99.    A string after ">" is the name of the following fragment.
  100.    Only the five upper-case letters A, C, G, T and N are allowed
  101.    to appear in fragment data. No other characters are allowed.
  102.    A common mistake is the use of lower case letters in a fragment.
  103.  
  104.    To run the program, type a command of form
  105.  
  106.     cap  file_of_fragments  
  107.  
  108.    The output goes to the terminal screen. So redirection of the
  109.    output into a file is necessary. The output consists of three parts:
  110.    overview of contigs at fragment level, detailed display of contigs
  111.    at nucleotide level, and consensus sequences.
  112.    '+' = direct orientation; '-' = reverse complement
  113.    The output of CAP on the sample input data looks like:
  114.  
  115. #Contig 1
  116.  
  117. #G022uabh+(0)
  118. TATTTTAGAGACCCAAGTTTTTGACCTTTTCCATGTTTACATCAATCCTGTAGGTGATTG
  119. GGCAGCCATTTAAGTATTATTATAGACATTTTCACTATCCCATTAAAACCCTTTATGCCC
  120. ATACATCATAACACTACTTCCTACCCATAAGCTCCTTTTAACTTGTTAAAGTCTTGCTTG
  121. AATTAAAGACTTGTTTAAACACAAAA-TTTAGACTTTTACTCAACAAAAGTGATTGATTG
  122. ATTGATTGATTGATTGAT
  123. #G028uaah+(145)
  124. CATAAGCTCCTTTTAACTTGTTAAAGTCTTGCTTGAATTAAAGACTTGTTTAAACACAAA
  125. A-TTTAGACTTTTACTCAACAAAAGTGATTGATTGATTGATTGATTGATTGATGGTTTAC
  126. AGTAGGACTTCATTCTAGTCATTATAGCTGCTGGCAGTATAACTGGCCAGCCTTTAATAC
  127. ATTGCTGCTTAGAGTCAAAGCATGTACTTAGAGTTGGTATGATTTATCTTTTTGGTCTTC
  128. TATAGCCTCCTTCCCCATCCC-ATCAGTCT
  129. #G019uabh+(120)
  130. ATACATCATAACACTACTTCCTACCCATAAGCTCCTTTTAACTTGTTAAAGTCTTGCTTG
  131. AATTAAAGACTTGTTTAAACACAAAAATTTAGAGTTTTACTCAACAAAAGTGATTGATTG
  132. ATTGATTGATTGATTGATGGTTTACAGTAGGACTTCATTCTAGTCATTATAGCTGCTGGC
  133. AGTATAACTGGCCAGCCTTTAATACATTGCTGCTTAGAGTCAAAGCATGTACTTAGAGTT
  134. GGTATGATTTATCTTTTTGGTCTTCTATAGCCTCCTTCCCCATCCCCATCAGTCTTAATC
  135. AGTCTTGTTACGTTATGACT-AATCTTTGGGGATTGTGCAGAATGTTATTTTAGATAAGC
  136. AAAA-CGAGCAAAAT-GGGGAGTT-A-CTT-A-TATTT-CTTT-AAA--GC
  137. #G023uabh-(426)
  138. GTAACGT-ATGA-TCAATCTTTGGGGATTGTGCAGAATGT-ATTTTAGATAAGCAAAAAC
  139. GAGCAAAATAGGGGAGTTTAACTTTAATATTTTCTTTTAAAAAGCATTTCATGTTATAAG
  140. ATCAATTCTGAGTGGTAGAAATGCTTTGACATTTTATTTCCATTTTCTACTTTTAGTTTT
  141. TTTCCTATTTGTTTAAGATCGTAGAGGATTATTAAGCTGAACTCCTCAACTGATAAAAAG
  142. CATGACATCTTAAACATAAGCAAAGCATATTTTTAGGTTAATTTTCACATAGAAAACAGT
  143. TTATTTTATGTGAAATTCTATGTAGATATACTATTTTTTTGGTATTTATT
  144. #G006uaah-(496)
  145. GGGGAGTTTAACTTTAATATTTTCTTTTAAAAAGCATTTCATGTTATAAGATCAATTCTG
  146. AGTGGTAGAAATGCTTTGACATTTTATTTCCATTTTCTACTTTTAGTTTTTTTCCTATTT
  147. GTTTAAGATCTTAGAGGATTATTAAGCTGAACTCCTCAACTGATAAAAAGCATGACATCT
  148. TAAACATAAGCAAAGCATATNNT-AGGTTAATTTTCACATAGAAAACAGTTTATTTTATG
  149. T
  150.  
  151.  
  152.  
  153. Slight modifications by S. Smith on Mon Feb 17 10:18:34 EST 1992.
  154. These changes allow for command line arguements for several
  155. of the hard coded parameters, as well as a slight modification to
  156. the output routine to support GDE format.  Changes are commented
  157. as:
  158.  
  159. Mod by S.S.
  160.  
  161. */
  162.  
  163.  
  164. #pragma segment Cap2
  165.  
  166. #include   <stdio.h>
  167. #include   <stdlib.h>
  168.  
  169.  
  170. int OVERLEN;        /* Minimum length of any overlap */
  171. float PERCENT;      /* Minimum identity percentage of any overlap */
  172. #define  CUTOFF    50            /* cutoff score for overlap or containment */
  173. #define  DELTA     9.0            /* Jump increment in check for overlap */
  174. #define  MATCH     2            /* score of a match */
  175. #define  MISMAT   -6            /* score of a mismatch */
  176. #define  LTMISM   -3            /* light score of a mismatch */
  177. #define  OPEN      0            /* gap open penalty */
  178. #define  EXTEND    4            /* gap extension penalty */
  179. #define  LTEXTEN   2            /* light gap extension penalty */
  180. #define  POS5      30            /* Sequencing errors often occur before base POS5 */
  181. #define  POS3      475            /* Sequencing errors often occur after base POS3 */
  182. #define  LINELEN   60        /* length of one printed line */
  183. #define  NAMELEN   20        /* length of printed fragment name */
  184. #define  TUPLELEN  4            /* Maximum length of words for lookup table */
  185. #define  SEQLEN    2000            /* initial size of an array for an output fragment */
  186.  
  187. static int over_len;
  188. static float per_cent;
  189. typedef struct OVERLAP        /* of 5' and 3' segments */
  190. {  
  191.     int    number;    /* index of 3' segment   */
  192.     int    host;        /* index of 5' segment   */
  193.     int    ind;        /* used in reassembly */
  194.     int    stari;    /* start position of 5' suffix */
  195.     int    endi;        /* end position of 5' suffix   */
  196.     int    starj;    /* start position of 3' prefix */
  197.     int    endj;        /* end position of 3' prefix   */
  198.     short   orienti;    /* orientation of 5' segment: 0= rev. */
  199.     short   orientj;    /* orientation of 3' segment: 0= rev. */
  200.     int    score;    /* score of overlap alignment */
  201.     int    length;    /* length of alignment      */
  202.     int    match;    /* number of matches in alignment */
  203.     short   kind;    /* 0 = containment; 1 = overlap   */
  204.     int    *script;    /* script for constructing alignment */
  205.     struct  OVERLAP  *next; 
  206. } over, *overptr;
  207. struct SEG
  208. {  
  209.     char    *name;    /* name string */
  210.     int     len;        /* length of segment name */
  211.     char    *seq;    /* segment sequence */
  212.     char    *rev;    /* reverse complement */
  213.     int     length;    /* length of sequence */
  214.     short   kind;    /* 0 = contain; 1 = overlap */
  215.     int     *lookup;    /* lookup table */
  216.     int     *pos;    /* location list */
  217.     overptr list;    /* list of overlapping edges */
  218. } *segment;    /* array of segment records */
  219. int   seg_num;                  /* The number of segments   */
  220. overptr   *edge;        /* set of overlapping edges */
  221. int   edge_num;                 /* The number of overlaps */
  222. struct CONS            /* 1 = itself; 0 = reverse complement */
  223. {  
  224.     short   isfive[2];    /* is 5' end free? */
  225.     short   isthree[2];    /* is 3' end free? */
  226.     short   orient[2];    /* orientation of 3' segment */
  227.     int     group;    /* contig serial number */
  228.     int     next[2];    /* pointer to 3' adjacent segment */
  229.     int     other;    /* the other end of the contig   */
  230.     int     child;    /* for the containment tree   */
  231.     int     brother;
  232.     int     father;
  233.     overptr node[2];    /* pointers to overlapping edges */
  234. } *contigs;  /* list of contigs */
  235. struct TTREE            /* multiple alignment tree */
  236. {  
  237.     short   head;    /* head = 1 means the head of a contig */
  238.     short   orient;    /* orientation */
  239.     int     begin;    /* start position of previous segment */
  240.     int     *script;    /* alignment script */
  241.     int     size;    /* length of script */
  242.     int     next;    /* list of overlap segments */
  243.     int     child;    /* list of child segments */
  244.     int     brother;    /* list of sibling segments */
  245. }  *mtree;
  246. int   vert[128];        /* converted digits for 'A','C','G','T' */
  247. int   vertc[128];        /* for reverse complement */
  248. int   tuple;            /* tuple = TUPLELEN - 1 */
  249. int  base[TUPLELEN];        /* locations of a lookup table */
  250. int  power[TUPLELEN];        /* powers of 4 */
  251. typedef struct OUT
  252. { 
  253.     char  *line;        /* one print line */
  254.     char  *a;        /* pointer to slot in line */
  255.     char  c;        /* current char */
  256.     char  *seq;        /* pointer to sequence */
  257.     int   length;        /* length of segment */
  258.     int   id;        /* index of segment */
  259.     int   *s;        /* pointer to script vector */
  260.     int   size;        /* size of script vector */
  261.     int   op;        /* current operation */
  262.     char  name[NAMELEN+2];/* name of segment */
  263.     short done;        /* indicates if segment is done */
  264.     int   loc;        /* position of next segment symbol */
  265.     char  kind;        /* type of next symbol of segment */
  266.     char  up;        /* type of upper symbol of operation */
  267.     char  dw;        /* type of lower symbol of operation */
  268.     int   offset;        /* relative to consensus sequence */
  269.     int   linesize;    /* size of array line */
  270.     struct  OUT *child;    /* pointer to child subtree */
  271.     struct  OUT *brother;    /* pointer to brother node */
  272.     struct  OUT *link;    /* for operation linked list */
  273.     struct  OUT *father;    /* pointer to father node */
  274. }  row, *rowptr;    /* node for segment */
  275. rowptr *work;            /* a set of working segments */
  276. rowptr head, tail;        /* first and last nodes of op list */
  277. struct VX
  278. {  
  279.     int     id;            /* Segment index */
  280.     short   kind;        /* overlap or containment */
  281.     overptr list;        /* list of overlapping edges */
  282. } *piece;        /* array of segment records */
  283. char  *allconsen, *allconpt;        /* Storing consensus sequences */
  284.  
  285. /* dgg patches for mac use */
  286.  
  287. char
  288. nuccomp( symbol)  char symbol;
  289. {
  290.     switch ( symbol ) { 
  291.         case 'A' : return  'T'; 
  292.         case 'a' : return  't'; 
  293.         case 'T' : return  'A'; 
  294.         case 't' : return  'a'; 
  295.         case 'C' : return  'G'; 
  296.         case 'c' : return  'g'; 
  297.         case 'G' : return  'C'; 
  298.         case 'g' : return  'c'; 
  299.         default  : return  symbol; 
  300.         }
  301. }
  302.  
  303. char *ckalloc(int amount); /* for borland c */
  304.  
  305.  
  306. typedef    char    cvec[128];
  307.  
  308. FILE* fout;    /* output file */
  309.  
  310.  
  311. #ifdef MACINTOSH
  312. // use with ChildApp.c
  313.  
  314. int RealMain( int argc, char *argv[])
  315.  
  316. #else
  317.  
  318. int main(argc, argv) 
  319. int argc; char *argv[];
  320. #endif
  321. { 
  322.     int   M;                    /* Sequence length */
  323. /* dgg -- make smaller arrays for mac */
  324. /*    int   V[128][128], Q,R;         
  325.     int   V1[128][128], R1;     
  326. */
  327.     cvec  *V, *V1;
  328.     short    Q,R, R1;     
  329.     int   total;                    /* Total of segment lengths */
  330.     int   number;                         /* The number of segments   */
  331.     char  *sequence;            /* Storing all segments     */
  332.     char  *reverse;            /* Storing all reverse segments */
  333.     int   symbol, prev;            /* The next character         */
  334.     FILE *Ap, *ckopen();                  /* For the sequence file      */
  335.     char *ckalloc();            /* space-allocating function  */
  336.     register int  i, j, k;        /* index variables          */
  337.     /* Mod by S.S. */
  338.     int jj;
  339.     short  heading;            /* 1: heading; 0: sequence    */
  340.     
  341.     /*
  342. *    Mod by S.S. & dgg
  343. *
  344.     if ( argc != 2 )
  345.         fatalf("The proper form of command is: \n%s file_of_fragments",
  346.              argv[0]);
  347. */
  348.     if ( argc < 5 ) {
  349.         short i;
  350.         for (i= 0; i<argc; i++)  fprintf(stderr, "arg[%d]=%s \n", i, argv[i]);
  351.         fatalf("usage: %s file_of_fragments output_file MIN_OVERLAP PERCENT_MATCH",
  352.             argv[0]);
  353.          }
  354.     sscanf(argv[3],"%d",&OVERLEN);
  355.     sscanf(argv[4],"%d",&jj);
  356.     PERCENT = (float)jj/100.0;
  357.     if(PERCENT < 0.25) PERCENT = 0.25;
  358.     if(PERCENT > 1.0) PERCENT = 1.0;
  359.     if(OVERLEN < 1) OVERLEN = 1;
  360.     if(OVERLEN > 100) OVERLEN = 100;
  361.  
  362.     if (argv[2] == "stdout" || argv[2] == "-") fout= NULL;
  363.     else fout = fopen( argv[2], "w");
  364.     if (!fout) fout = stdout;
  365.  
  366.     /* determine number of segments and total lengths */
  367.     j = 0;
  368.  
  369.     Ap = ckopen(argv[1], "r");
  370.     prev = '\n';
  371.     for (total = 3, number = 0; ( symbol = fgetc(Ap)) != EOF ; total++ )
  372.     { 
  373.         if ( symbol == '>' && prev == '\n' )
  374.             number++;
  375.         prev = symbol;
  376.     }
  377.     (void)    fclose(Ap);
  378.  
  379.     /* dgg - allocate space for arrays */
  380.     /* for (i=0; i<128; i++) V[i]= ( char *) ckalloc( 128 * sizeof(char)); NO NO */
  381.     V= ( cvec *) ckalloc( 128 * 128 * sizeof(char)); /* !! yes */
  382.     V1= ( cvec *) ckalloc( 128 * 128 * sizeof(char));
  383.  
  384.     total += number * 20;
  385.     /* allocate space for segments */
  386.     sequence = ( char * ) ckalloc( total * sizeof(char));
  387.     reverse = ( char * ) ckalloc( total * sizeof(char));
  388.     allconpt = allconsen = ( char * ) ckalloc( total * sizeof(char));
  389.     segment = ( struct SEG * ) ckalloc( number * sizeof(struct SEG));
  390.  
  391.     /* read the segments into sequence */
  392.     M = 0;
  393.     Ap = ckopen(argv[1], "r");
  394.     number = -1;
  395.     heading = 0;
  396.     prev = '\n';
  397.     for ( i = 0, k = total ; ( symbol = fgetc(Ap)) != EOF ; )
  398.     { 
  399.         if ( symbol != '\n' )
  400.         { 
  401.             sequence[++i] = symbol;
  402.             /* reverse[--k] = nuccomp( symbol); */
  403.             switch ( symbol )
  404.             { 
  405.             case 'A' : 
  406.                 reverse[--k] = 'T'; 
  407.                 break;
  408.             case 'a' : 
  409.                 reverse[--k] = 't'; 
  410.                 break;
  411.             case 'T' : 
  412.                 reverse[--k] = 'A'; 
  413.                 break;
  414.             case 't' : 
  415.                 reverse[--k] = 'a'; 
  416.                 break;
  417.             case 'C' : 
  418.                 reverse[--k] = 'G'; 
  419.                 break;
  420.             case 'c' : 
  421.                 reverse[--k] = 'g'; 
  422.                 break;
  423.             case 'G' : 
  424.                 reverse[--k] = 'C'; 
  425.                 break;
  426.             case 'g' : 
  427.                 reverse[--k] = 'c'; 
  428.                 break;
  429.             default  : 
  430.                 reverse[--k] = symbol; 
  431.                 break;
  432.             }
  433.         }
  434.         if ( symbol == '>' && prev == '\n' )
  435.         { 
  436.             heading = 1;
  437.             if ( number >= 0 )
  438.             { 
  439.                 segment[number].length = i - j - 1;
  440.                 segment[number].rev = &(reverse[k]);
  441.                 if ( i - j - 1 > M ) M = i - j -1;
  442.             }
  443.             number++;
  444.             j = i;
  445.             segment[number].name = &(sequence[i]);
  446.             segment[number].kind = 1;
  447.             segment[number].list = 0;
  448.         }
  449.         if ( heading && symbol == '\n' )
  450.         { 
  451.             heading = 0;
  452.             segment[number].len = i - j;
  453.             segment[number].seq = &(sequence[i]);
  454.             j = i;
  455.         }
  456.  
  457.         prev = symbol;
  458.     }
  459.     segment[number].length = i - j;
  460.     reverse[--k] = '>';
  461.     segment[number].rev = &(reverse[k]);
  462.     if ( i - j > M ) M = i - j;
  463.     seg_num = ++number;
  464.     (void)    fclose(Ap);
  465.  
  466.     Q = OPEN;
  467.     R = EXTEND;
  468.     R1 = LTEXTEN;
  469.     /* set match and mismatch weights */
  470.     for ( i = 0; i < 128 ; i++ )
  471.         for ( j = 0; j < 128 ; j++ )
  472.             if ((i|32) == (j|32) )
  473.                 V[i][j] = V1[i][j] = MATCH;
  474.             else
  475.             { 
  476.                 V[i][j] = MISMAT;
  477.                 V1[i][j] = LTMISM;
  478.             }
  479.     for ( i = 0; i < 128 ; i++ )
  480.         V['N'][i] = V[i]['N'] = MISMAT + 1;
  481.     V1['N']['N'] = MISMAT + 1;
  482.  
  483.     over_len = OVERLEN;
  484.     per_cent = PERCENT;
  485.     edge_num = 0;
  486.     fprintf(stderr,"\nINIT"); INIT(M);
  487.     fprintf(stderr,"\nMAKE"); MAKE();
  488.     fprintf(stderr,"\nPAIR"); PAIR(V,V1,Q,R,R1);
  489.     fprintf(stderr,"\nASSEM"); ASSEM();
  490.     fprintf(stderr,"\nREPAIR"); REPAIR();
  491.     fprintf(stderr,"\nFORM_TREE"); FORM_TREE();
  492.     /* GRAPH(); */
  493.     fprintf(stderr,"\nSHOW\n"); SHOW();  
  494.     fclose(fout); /* dgg */
  495. }
  496.  
  497. static char  (*varray)[128];            /* substitution scores */
  498. static int  q, r;            /* gap penalties */
  499. static int  qr;                /* qr = q + r */
  500. static char  (*v1)[128];            /* light substitution scores */
  501. static int  r1;                /* light extension penalty */
  502. static int  qr1;            /* qr1 = q + r1 */
  503.  
  504. static int   SCORE;
  505. static int   STARI;
  506. static int   STARJ;
  507. static int   ENDI;
  508. static int   ENDJ;
  509.  
  510. static int  *CC, *DD;            /* saving matrix scores */
  511. static int  *RR, *SS;             /* saving start-points */
  512. static int   *S;            /* saving operations for diff */
  513.  
  514. /* The following definitions are for function diff() */
  515.  
  516. int   diff();
  517. static int   zero = 0;                /* int  type zero        */
  518.  
  519. #define gap(k)  ((k) <= 0 ? 0 : q+r*(k))    /* k-symbol indel score */
  520.  
  521. static int  *sapp;                /* Current script append ptr */
  522. static int   last;                /* Last script op appended */
  523.  
  524. static int  no_mat;                 /* number of matches */
  525.  
  526. static int  no_mis;                 /* number of mismatches */
  527.  
  528. static int  al_len;                 /* length of alignment */
  529. /* Append "Delete k" op */
  530. #define DEL(k)                \
  531. { al_len += k;                \
  532.   if (last < 0)                \
  533.     last = sapp[-1] -= (k);        \
  534.   else                    \
  535.     last = *sapp++ = -(k);        \
  536. }
  537. /* Append "Insert k" op */
  538. #define INS(k)                \
  539. { al_len += k;                \
  540.   if (last < 0)                \
  541.     { sapp[-1] = (k); *sapp++ = last; }    \
  542.   else                    \
  543.     last = *sapp++ = (k);        \
  544. }
  545.  
  546. /* Append "Replace" op */
  547. #define REP                 \
  548. { last = *sapp++ = 0;             \
  549.   al_len += 1;                \
  550. }
  551.  
  552. INIT(M) int  M;
  553. { 
  554.     register  int   j;            /* row and column indices */
  555.     char *ckalloc();            /* space-allocating function */
  556.  
  557.     /* allocate space for all vectors */
  558.     j = (M + 1) * sizeof(int );
  559.     CC = ( int  * ) ckalloc(j);
  560.     DD = ( int  * ) ckalloc(j);
  561.     RR = ( int  * ) ckalloc(j);
  562.     SS = ( int  * ) ckalloc(j);
  563.     S = ( int  * ) ckalloc(2 * j);
  564. }
  565.  
  566. /* Make a lookup table for words of lengths up to TUPLELEN in each sequence.
  567.    The value of a word is used as an index to the table.  */
  568. MAKE()
  569. { 
  570.     int   hash[TUPLELEN];        /* values of words of lengths up to TUPLELEN */
  571.     int   *table;            /* pointer to a lookup table */
  572.     int   *loc;            /* pointer to a table of sequence locations */
  573.     int   size;            /* size of a lookup table */
  574.     int   limit, digit, step;    /* temporary variables  */
  575.     register  int  i, j, k, p, q;    /* index varibles */
  576.     char *ckalloc();        /* space-allocating function */
  577.     char *A;            /* pointer to a sequence */
  578.     int  M;            /* length of a sequence */
  579.  
  580.     tuple = TUPLELEN - 1;
  581.     for ( i = 0; i < 128; i++ )
  582.         vert[i] = 4;
  583.     vert['A'] = vert['a'] = 0;
  584.     vert['C'] = vert['c'] = 1;
  585.     vert['G'] = vert['g'] = 2;
  586.     vert['T'] = vert['t'] = 3;
  587.     vertc['A'] = vertc['a'] = 3;
  588.     vertc['C'] = vertc['c'] = 2;
  589.     vertc['G'] = vertc['g'] = 1;
  590.     vertc['T'] = vertc['t'] = 0;
  591.     for ( i = 0, j = 1, size = 0; i <= tuple ; i++, j *= 4 )
  592.     {  
  593.         base[i] = size;
  594.         power[i] = j;
  595.         size = ( size + 1 ) * 4;
  596.     }
  597.     for ( j = 0; j <= tuple; j++ )
  598.         hash[j] = 0;
  599.     /* make a lookup table for each sequence */
  600.     for ( i = 0; i < seg_num ; i++ )
  601.     { 
  602.         A = segment[i].seq;
  603.         M = segment[i].length;
  604.         table = segment[i].lookup = (int  * ) ckalloc(size * sizeof(int ));
  605.         loc = segment[i].pos = (int  * ) ckalloc((M + 1) * sizeof(int ));
  606.         for ( j = 0; j < size; j++ )
  607.             table[j] = 0;
  608.         for ( k = 0, j = 1; j <= M; j++ )
  609.             if ( ( digit = vert[A[j]] ) != 4 )
  610.             { 
  611.                 for ( p = tuple; p > 0; p-- )
  612.                     hash[p] = 4 * (hash[p-1] + 1) + digit;
  613.                 hash[0] = digit;
  614.                 step = j - k;
  615.                 limit = tuple <= step ? tuple : step;
  616.                 for ( p = 0; p < limit; p++ )
  617.                     if ( ! table[(q = hash[p])] )  table[q] = 1;
  618.                 if ( step > tuple )
  619.                 { 
  620.                     loc[(p = j - tuple)] = table[(q = hash[tuple])];
  621.                     table[q] = p;
  622.                 }
  623.             }
  624.             else
  625.                 k = j;
  626.     }
  627. }
  628.  
  629. /* Perform pair-wise comparisons of sequences. The sequences not
  630.    satisfying the necessary condition for overlap are rejected quickly.
  631.    Those that do satisfy the condition are verified with a dynamic
  632.    programming algorithm to see if overlaps exist. */
  633. PAIR(V,V1,Q,R,R1) char  V[][128],V1[][128],Q,R,R1;
  634. { 
  635.     int  endi, endj, stari, starj;    /* endpoint and startpoint */
  636.  
  637.     short orienti, orientj;        /* orientation of segments */
  638.     short iscon;                /* containment condition   */
  639.     int   score;               /* the max score */
  640.     int  count, limit;            /* temporary variables     */
  641.  
  642.     register  int   i, j, d;        /* row and column indices  */
  643.     char *ckalloc();            /* space-allocating function */
  644.     int  rl, cl;
  645.     char *A, *B;
  646.     int  M, N;
  647.     overptr node1;
  648.     int  total;                /* total number of pairs */
  649.     int  hit;                /* number of pairs satisfying cond. */
  650.     int  CHECK();
  651.  
  652.     varray = V;
  653.     q = Q;
  654.     r = R;
  655.     qr = q + r;
  656.     v1 = V1;
  657.     r1 = R1;
  658.     qr1 = q + r1;
  659.     total = hit = 0;
  660.     limit = 2 * ( seg_num - 1 );
  661.     for ( orienti = 0, d = 0; d < limit ; d++ )
  662.     { 
  663.         i = d / 2;
  664.         orienti = 1 - orienti;
  665.         A = orienti ? segment[i].seq : segment[i].rev;
  666.         M = segment[i].length;
  667.         for ( j = i+1; j < seg_num ; j++ )
  668.         { 
  669.             B = segment[j].seq;
  670.             orientj = 1;
  671.             N = segment[j].length;
  672.             total += 1;
  673.             if ( CHECK(orienti, i, j) )
  674.             { 
  675.                 hit += 1;
  676.                 SCORE = 0;
  677.                 big_pass(A,B,M,N,orienti,orientj);
  678.                 if ( SCORE )
  679.                 { 
  680.                     score = SCORE;
  681.                     stari = ++STARI;
  682.                     starj = ++STARJ;
  683.                     endi = ENDI;
  684.                     endj = ENDJ;
  685.                     rl = endi - stari + 1;
  686.                     cl = endj - starj + 1;
  687.                     sapp = S;
  688.                     last = 0;
  689.                     al_len = no_mat = no_mis = 0;
  690.                     (void) diff(&A[stari]-1, &B[starj]-1,rl,cl,q,q);
  691.                     iscon = stari == 1 && endi == M || starj == 1 && endj == N;
  692.                     if ( no_mat >= al_len * per_cent &&
  693.                         ( al_len >= over_len || iscon ) )
  694.                     { 
  695.                         node1 = ( overptr ) ckalloc( (int ) sizeof(over));
  696.                         if ( iscon )
  697.                             node1->kind = 0;        /* containment */
  698.                         else
  699.                         { 
  700.                             node1->kind = 1; 
  701.                             edge_num++; 
  702.                         }    /* overlap */
  703.                         if ( endi == M && ( endj != N || starj != 1 ) ) /*i is 5'*/
  704.                         { 
  705.                             node1->number = j;
  706.                             node1->host = i;
  707.                             node1->stari = stari;
  708.                             node1->endi = endi;
  709.                             node1->orienti = orienti;
  710.                             node1->starj = starj;
  711.                             node1->endj = endj;
  712.                             node1->orientj = orientj;
  713.                         }
  714.                         else    /* j is 5' */
  715.                         { 
  716.                             node1->number = i;
  717.                             node1->host = j;
  718.                             node1->stari = starj;
  719.                             node1->endi = endj;
  720.                             node1->orienti = orientj;
  721.                             node1->starj = stari;
  722.                             node1->endj = endi;
  723.                             node1->orientj = orienti;
  724.                         }
  725.                         node1->score = score;
  726.                         node1->length = al_len;
  727.                         node1->match = no_mat;
  728.                         count = node1->number == i ? j : i;
  729.                         node1->next = segment[count].list;
  730.                         segment[count].list = node1;
  731.                         if ( ! node1->kind )
  732.                             segment[count].kind = 0;
  733.                     }
  734.                 }
  735.             }
  736.         }
  737.     }
  738. }
  739.  
  740. /* Return 1 if two sequences satisfy the necessary condition for overlap,
  741.    and 0 otherwise. Parameters first and second are the indices of segments,
  742.    and parameter orient indicates the orientation of segment first.  */
  743. int  CHECK(orient,first,second) short orient; 
  744. int  first, second;
  745. { 
  746.     int   limit, bound;        /* maximum number of jumps */
  747.     int   small;            /* smaller of limit and bound */
  748.     float  delta;            /* cutoff factor   */
  749.     float cut;            /* cutoff score */
  750.     register  int   i;        /* index variable */
  751.     int  t, q;            /* temporary variables */
  752.     int  JUMP();
  753.     int  RJUMP();
  754.     int  JUMPC();
  755.     int  RJUMPC();
  756.  
  757.     delta = DELTA;
  758.     if ( orient )
  759.         limit = JUMP(CC, first, second, 1);
  760.     else
  761.         limit = JUMPC(CC, first, second);
  762.     bound = RJUMP(DD, second, first, orient);
  763.     small = limit <= bound ? limit : bound;
  764.     cut = 0;
  765.     for ( i = 1; i <= small; i++ )
  766.     { 
  767.         if ( (t = DD[i] - 1) >= over_len && t >= cut
  768.             && (q = CC[i] - 1) >= over_len && q >= cut )
  769.             return (1);
  770.         cut += delta;
  771.     }
  772.     if (DD[bound] >= delta*(bound-1) || CC[limit] >= delta*(limit-1))
  773.         return (1);
  774.     limit = JUMP(CC, second, first, orient);
  775.     if ( orient )
  776.         bound = RJUMP(DD, first, second, 1);
  777.     else
  778.         bound = RJUMPC(DD, first, second);
  779.     small = limit <= bound ? limit : bound;
  780.     cut = 0;
  781.     for ( i = 1; i <= small; i++ )
  782.     { 
  783.         if ( (t = DD[i] - 1) >= over_len && t >= cut
  784.             && (q = CC[i] - 1) >= over_len && q >= cut )
  785.             return (1);
  786.         cut += delta;
  787.     }
  788.     return (0);
  789. }
  790.  
  791. /* Compute a vector of lengths of jumps */
  792. int  JUMP(H,one,two,orient) int  H[], one, two; 
  793. short orient;
  794. { 
  795.     char *A, *B;            /* pointers to two sequences */
  796.     int  M, N;            /* lengths of two sequences */
  797.     int  *table;            /* pointer to a lookup table */
  798.     int  *loc;            /* pointer to a location table */
  799.     int  value;            /* value of a word */
  800.     int  maxd;            /* maximum length of an identical diagonal */
  801.     int  d;            /* length of current identical diagonal */
  802.     int  s;            /* length of jumps */
  803.     int  k;            /* number of jumps */
  804.  
  805.     register int  i, j, p;    /* index variables */
  806.  
  807.     A = segment[one].seq;
  808.     M = segment[one].length;
  809.     table = segment[one].lookup;
  810.     loc = segment[one].pos;
  811.     B = orient ? segment[two].seq : segment[two].rev;
  812.     N = segment[two].length;
  813.     for ( s = 1, k = 1; s <= N ; k++ )
  814.     { 
  815.         maxd = 0;
  816.         for ( value = -1, d = 0, j = s; d <= tuple && j <= N; j++, d++)
  817.         { 
  818.             if ( vert[B[j]] == 4 ) break;
  819.             value = 4 * (value + 1) + vert[B[j]];
  820.             if ( ! table[value] ) break;
  821.         }
  822.         if ( d > tuple )
  823.         { 
  824.             for ( p = table[value]; p ; p = loc[p] )
  825.             { 
  826.                 d = tuple + 1;
  827.                 for ( i = p+d, j = s+d; i <= M && j <= N; i++, j++, d++ )
  828.                     if ( A[i] != B[j] && vert[A[i]] != 4 && vert[B[j]] != 4 ) break;
  829.                 if ( maxd < d )
  830.                     maxd = d;
  831.                 if ( j > N ) break;
  832.             }
  833.         }
  834.         else 
  835.             maxd = d;
  836.         s += maxd + 1;
  837.         H[k] = s;
  838.     }
  839.     return (k - 1);
  840. }
  841.  
  842. /* Compute a vector of lengths of jumps for reverse complement of one */
  843. int  JUMPC(H,one,two) int  H[], one, two;
  844. { 
  845.     char *A, *B;            /* pointers to two sequences */
  846.     int  M, N;            /* lengths of two sequences */
  847.     int  *table;            /* pointer to a lookup table */
  848.     int  *loc;        /* pointer to a location table */
  849.     int  value;            /* value of a word */
  850.     int  maxd;            /* maximum length of an identical diagonal */
  851.     int  d;            /* length of current identical diagonal */
  852.     int  s;            /* length of jumps */
  853.     int  k;            /* number of jumps */
  854.  
  855.     register  int   i, j, p;    /* index variables */
  856.  
  857.     A = segment[one].rev;
  858.     M = segment[one].length;
  859.     table = segment[one].lookup;
  860.     loc = segment[one].pos;
  861.     B = segment[two].seq;
  862.     N = segment[two].length;
  863.     for ( s = 1, k = 1; s <= N ; k++ )
  864.     { 
  865.         maxd = 0;
  866.         for ( value = 0, d = 0, j = s; d <= tuple && j <= N; j++, d++)
  867.         { 
  868.             if ( vert[B[j]] == 4 ) break;
  869.             value += vertc[B[j]] * power[d];
  870.             if ( ! table[value+base[d]] ) break;
  871.         }
  872.         if ( d > tuple )
  873.         { 
  874.             for ( p = table[value+base[tuple]]; p ; p = loc[p] )
  875.             { 
  876.                 d = tuple + 1;
  877.                 for ( i = M+2-p, j = s+d; i <= M && j <= N; i++, j++, d++ )
  878.                     if ( A[i] != B[j] && vert[A[i]] != 4 && vert[B[j]] != 4 ) break;
  879.                 if ( maxd < d )
  880.                     maxd = d;
  881.                 if ( j > N ) break;
  882.             }
  883.         }
  884.         else 
  885.             maxd = d;
  886.         s += maxd + 1;
  887.         H[k] = s;
  888.     }
  889.     return (k - 1);
  890. }
  891.  
  892. /* Compute a vector of lengths of reverse jumps */
  893. int  RJUMP(H,one,two,orient) int  H[], one, two; 
  894. short orient;
  895. { 
  896.     char *A, *B;            /* pointers to two sequences */
  897.     int  N;            /* length of a sequence */
  898.     int  *table;            /* pointer to a lookup table */
  899.     int  *loc;        /* pointer to a location table */
  900.     int  value;            /* value of a word */
  901.     int  maxd;            /* maximum length of an identical diagonal */
  902.     int  d;            /* length of current identical diagonal */
  903.     int  s;            /* length of jumps */
  904.     int  k;            /* number of jumps */
  905.  
  906.     register  int   i, j, p;    /* index variables */
  907.  
  908.     A = segment[one].seq;
  909.     table = segment[one].lookup;
  910.     loc = segment[one].pos;
  911.     B = orient ? segment[two].seq : segment[two].rev;
  912.     N = segment[two].length;
  913.     for ( s = 1, k = 1; s <= N ; k++ )
  914.     { 
  915.         maxd = 0;
  916.         for (value = 0, d = 0, j = N+1-s; d <= tuple && j >= 1; j--, d++)
  917.         { 
  918.             if ( vert[B[j]] == 4 ) break;
  919.             value += vert[B[j]] * power[d];
  920.             if ( ! table[value+base[d]] ) break;
  921.         }
  922.         if ( d > tuple )
  923.         { 
  924.             for ( p = table[value+base[tuple]]; p ; p = loc[p] )
  925.             { 
  926.                 d = tuple + 1;
  927.                 for ( i = p-1, j = N+1-s-d; i >= 1 && j >= 1; i--, j--, d++ )
  928.                     if ( A[i] != B[j] && vert[A[i]] != 4 && vert[B[j]] != 4 ) break;
  929.                 if ( maxd < d )
  930.                     maxd = d;
  931.                 if ( j < 1 ) break;
  932.             }
  933.         }
  934.         else 
  935.             maxd = d;
  936.         s += maxd + 1;
  937.         H[k] = s;
  938.     }
  939.     return (k - 1);
  940. }
  941.  
  942. /* Compute a vector of lengths of reverse jumps for reverse complement */
  943. int  RJUMPC(H,one,two) int  H[], one, two;
  944. { 
  945.     char *A, *B;            /* pointers to two sequences */
  946.     int  M, N;            /* lengths of two sequences */
  947.     int  *table;            /* pointer to a lookup table */
  948.     int  *loc;        /* pointer to a location table */
  949.     int  value;            /* value of a word */
  950.     int  maxd;            /* maximum length of an identical diagonal */
  951.     int  d;            /* length of current identical diagonal */
  952.     int  s;            /* length of jumps */
  953.     int  k;            /* number of jumps */
  954.  
  955.     register  int   i, j, p;    /* index variables */
  956.  
  957.     A = segment[one].rev;
  958.     M = segment[one].length;
  959.     table = segment[one].lookup;
  960.     loc = segment[one].pos;
  961.     B = segment[two].seq;
  962.     N = segment[two].length;
  963.     for ( s = 1, k = 1; s <= N ; k++ )
  964.     { 
  965.         maxd = 0;
  966.         for (value = -1, d = 0, j = N+1-s; d <= tuple && j >= 1; j--, d++)
  967.         { 
  968.             if ( vert[B[j]] == 4 ) break;
  969.             value = 4 * (value + 1) + vertc[B[j]];
  970.             if ( ! table[value] ) break;
  971.         }
  972.         if ( d > tuple )
  973.         { 
  974.             for ( p = table[value]; p ; p = loc[p] )
  975.             { 
  976.                 d = tuple + 1;
  977.                 i = M - p - tuple;
  978.                 for ( j = N-s-tuple; i >= 1 && j >= 1; i--, j--, d++ )
  979.                     if ( A[i] != B[j] && vert[A[i]] != 4 && vert[B[j]] != 4 ) break;
  980.                 if ( maxd < d )
  981.                     maxd = d;
  982.                 if ( j < 1 ) break;
  983.             }
  984.         }
  985.         else 
  986.             maxd = d;
  987.         s += maxd + 1;
  988.         H[k] = s;
  989.     }
  990.     return (k - 1);
  991. }
  992.  
  993. /* Construct contigs */
  994. ASSEM()
  995. { 
  996.     char *ckalloc();        /* space-allocating function */
  997.     register  int   i, j, k;    /* index variables */
  998.     overptr  node1, x, y;        /* temporary pointer */
  999.     int   five, three;        /* indices of 5' and 3' segments */
  1000.     short   orienti;        /* orientation of 5' segment */
  1001.     short   orientj;        /* orientation of 3' segment */
  1002.     short   sorted;        /* boolean variable */
  1003.  
  1004.     contigs = ( struct CONS * ) ckalloc( seg_num * sizeof(struct CONS));
  1005.     for ( i = 0; i < seg_num; i++ )
  1006.     { 
  1007.         contigs[i].isfive[0] = contigs[i].isthree[0] = 1;
  1008.         contigs[i].isfive[1] = contigs[i].isthree[1] = 1;
  1009.         contigs[i].other = i;
  1010.         contigs[i].group = contigs[i].child = -1;
  1011.         contigs[i].brother = contigs[i].father = -1;
  1012.     }
  1013.     for ( i = 0; i < seg_num; i++ )
  1014.         if ( ! segment[i].kind )
  1015.             for ( ; ; )
  1016.             { 
  1017.                 for ( y = segment[i].list; y->kind; y = y->next )
  1018.                     ;
  1019.                 for ( x = y->next; x != 0; x = x->next )
  1020.                     if ( ! x->kind && x->score > y->score )
  1021.                         y = x;
  1022.                 for ( j = y->number; (k = contigs[j].father) != -1; j = k )
  1023.                     ;
  1024.                 if ( j != i )
  1025.                 { 
  1026.                     contigs[i].father = j = y->number;
  1027.                     contigs[i].brother = contigs[j].child;
  1028.                     contigs[j].child = i;
  1029.                     contigs[i].node[1] = y;
  1030.                     break;
  1031.                 }
  1032.                 else
  1033.                 { 
  1034.                     if ( segment[i].list->number == y->number )
  1035.                         segment[i].list = y->next;
  1036.                     else
  1037.                     { 
  1038.                         for ( x = segment[i].list; x->next->number != y->number ; )
  1039.                             x = x->next;
  1040.                         x->next = y->next;
  1041.                     }
  1042.                     for ( x = segment[i].list; x != 0 && x->kind; x = x->next )
  1043.                         ;
  1044.                     if ( x == 0 )
  1045.                     { 
  1046.                         segment[i].kind = 1;
  1047.                         break;
  1048.                     }
  1049.                 }
  1050.             }
  1051.     edge = ( overptr * ) ckalloc( edge_num * sizeof(overptr) );
  1052.     for ( j = 0, i = 0; i < seg_num; i++ )
  1053.         if ( segment[i].kind )
  1054.             for ( node1 = segment[i].list; node1 != 0; node1 = node1->next )
  1055.                 if ( segment[node1->number].kind )
  1056.                     edge[j++] = node1;
  1057.     edge_num = j;
  1058.     for ( i = edge_num - 1; i > 0; i-- )
  1059.     { 
  1060.         sorted = 1;
  1061.         for ( j = 0; j < i; j++ )
  1062.             if ( edge[j]->score < edge[j+1]->score )
  1063.             { 
  1064.                 node1 = edge[j];
  1065.                 edge[j] = edge[j+1];
  1066.                 edge[j+1] = node1;
  1067.                 sorted = 0;
  1068.             }
  1069.         if ( sorted )
  1070.             break;
  1071.     }
  1072.     for ( k = 0; k < edge_num; k++ )
  1073.     { 
  1074.         five = edge[k]->host;
  1075.         three = edge[k]->number;
  1076.         orienti = edge[k]->orienti;
  1077.         orientj = edge[k]->orientj;
  1078.         if ( contigs[five].isthree[orienti] &&
  1079.             contigs[three].isfive[orientj] && contigs[five].other != three )
  1080.         { 
  1081.             contigs[five].isthree[orienti] = 0;
  1082.             contigs[three].isfive[orientj] = 0;
  1083.             contigs[five].next[orienti] = three;
  1084.             contigs[five].orient[orienti] = orientj;
  1085.             contigs[five].node[orienti] = edge[k];
  1086.             contigs[three].isthree[(j = 1 - orientj)] = 0;
  1087.             contigs[five].isfive[(i = 1 - orienti)] = 0;
  1088.             contigs[three].next[j] = five;
  1089.             contigs[three].orient[j] = i;
  1090.             contigs[three].node[j] = edge[k];
  1091.             i = contigs[three].other;
  1092.             j = contigs[five].other;
  1093.             contigs[i].other = j;
  1094.             contigs[j].other = i;
  1095.         }
  1096.     }
  1097. }
  1098.  
  1099. REPAIR()
  1100. { 
  1101.     int  endi, endj, stari, starj;    /* endpoint and startpoint */
  1102.  
  1103.     short orienti, orientj;        /* orientation of segments */
  1104.     short isconi, isconj;            /* containment condition   */
  1105.     int   score;               /* the max score */
  1106.     int   i, j, f, d, e;            /* row and column indices  */
  1107.     char *ckalloc();            /* space-allocating function */
  1108.     char *A, *B;
  1109.     int  M, N;
  1110.     overptr node1;
  1111.     int   piece_num;                   /* The number of pieces */
  1112.     int   count, limit;
  1113.     int   number;
  1114.     int   hit;
  1115.  
  1116.     piece = ( struct VX * ) ckalloc( seg_num * sizeof(struct VX));
  1117.     for ( j = 0, i = 0; i < seg_num; i++ )
  1118.         if (segment[i].kind && (contigs[i].isfive[1] || contigs[i].isfive[0]))
  1119.             piece[j++].id = i;
  1120.     piece_num = j;
  1121.     for ( i = 0; i < piece_num; i++ )
  1122.     { 
  1123.         piece[i].kind = 1;
  1124.         piece[i].list = 0;
  1125.     }
  1126.     limit = 2 * ( piece_num - 1 );
  1127.     hit = number = 0;
  1128.     for ( orienti = 0, d = 0; d < limit ; d++ )
  1129.     { 
  1130.         i = piece[(e = d / 2)].id;
  1131.         orienti = 1 - orienti;
  1132.         A = orienti ? segment[i].seq : segment[i].rev;
  1133.         M = segment[i].length;
  1134.         for ( f = e+1; f < piece_num ; f++ )
  1135.         { 
  1136.             j = piece[f].id;
  1137.             B = segment[j].seq;
  1138.             orientj = 1;
  1139.             N = segment[j].length;
  1140.             SCORE = 0;
  1141.             hit++;
  1142.             big_pass(A,B,M,N,orienti,orientj);
  1143.             if ( SCORE > CUTOFF )
  1144.             { 
  1145.                 score = SCORE;
  1146.                 stari = ++STARI;
  1147.                 starj = ++STARJ;
  1148.                 endi = ENDI;
  1149.                 endj = ENDJ;
  1150.                 isconi = stari == 1 && endi == M;
  1151.                 isconj = starj == 1 && endj == N;
  1152.                 node1 = ( overptr ) ckalloc( (int ) sizeof(over));
  1153.                 if ( isconi || isconj )
  1154.                     node1->kind = 0;        /* containment */
  1155.                 else
  1156.                 { 
  1157.                     node1->kind = 1; 
  1158.                     number++; 
  1159.                 }    /* overlap */
  1160.                 if ( endi == M && ! isconj )     /*i is 5'*/
  1161.                 { 
  1162.                     node1->number = j;
  1163.                     node1->host = i;
  1164.                     node1->ind = f;
  1165.                     node1->stari = stari;
  1166.                     node1->endi = endi;
  1167.                     node1->orienti = orienti;
  1168.                     node1->starj = starj;
  1169.                     node1->endj = endj;
  1170.                     node1->orientj = orientj;
  1171.                 }
  1172.                 else    /* j is 5' */
  1173.                 { 
  1174.                     node1->number = i;
  1175.                     node1->host = j;
  1176.                     node1->ind = e;
  1177.                     node1->stari = starj;
  1178.                     node1->endi = endj;
  1179.                     node1->orienti = orientj;
  1180.                     node1->starj = stari;
  1181.                     node1->endj = endi;
  1182.                     node1->orientj = orienti;
  1183.                 }
  1184.                 node1->score = score;
  1185.                 count = node1->number == i ? f : e;
  1186.                 node1->next = piece[count].list;
  1187.                 piece[count].list = node1;
  1188.                 if ( ! node1->kind )
  1189.                     piece[count].kind = 0;
  1190.             }
  1191.         }
  1192.     }
  1193.     REASSEM(piece_num, number);
  1194. }
  1195.  
  1196. /* Construct contigs */
  1197. REASSEM(piece_num, number) int  piece_num, number;
  1198. { 
  1199.     char *ckalloc();        /* space-allocating function */
  1200.     int   i, j, k, d;        /* index variables */
  1201.     overptr  node1, x, y;        /* temporary pointer */
  1202.     int   five, three;        /* indices of 5' and 3' segments */
  1203.     short   orienti;        /* orientation of 5' segment */
  1204.     short   orientj;        /* orientation of 3' segment */
  1205.     short   sorted;        /* boolean variable */
  1206.  
  1207.     for ( d = 0; d < piece_num; d++ )
  1208.         if ( ! piece[d].kind )
  1209.             for ( i = piece[d].id ; ; )
  1210.             { 
  1211.                 for ( y = piece[d].list; y->kind; y = y->next )
  1212.                     ;
  1213.                 for ( x = y->next; x != 0; x = x->next )
  1214.                     if ( ! x->kind && x->score > y->score )
  1215.                         y = x;
  1216.                 for ( j = y->number; (k = contigs[j].father) != -1; j = k )
  1217.                     ;
  1218.                 if ( j != i && RECONCILE(y,&piece_num,&number) )
  1219.                 { 
  1220.                     contigs[i].father = j = y->number;
  1221.                     contigs[i].brother = contigs[j].child;
  1222.                     contigs[j].child = i;
  1223.                     contigs[i].node[1] = y;
  1224.                     segment[i].kind = 0;
  1225.                     break;
  1226.                 }
  1227.                 else
  1228.                 { 
  1229.                     if ( piece[d].list->number == y->number )
  1230.                         piece[d].list = y->next;
  1231.                     else
  1232.                     { 
  1233.                         for ( x = piece[d].list; x->next->number != y->number ; )
  1234.                             x = x->next;
  1235.                         x->next = y->next;
  1236.                     }
  1237.                     for ( x = piece[d].list; x != 0 && x->kind; x = x->next )
  1238.                         ;
  1239.                     if ( x == 0 )
  1240.                     { 
  1241.                         piece[d].kind = 1;
  1242.                         break;
  1243.                     }
  1244.                 }
  1245.             }
  1246.     if ( number > edge_num )
  1247.         edge = ( overptr * ) ckalloc( number * sizeof(overptr) );
  1248.     for ( j = 0, d = 0; d < piece_num; d++ )
  1249.         if ( piece[d].kind )
  1250.             for ( node1 = piece[d].list; node1 != 0; node1 = node1->next )
  1251.                 if ( piece[node1->ind].kind )
  1252.                     edge[j++] = node1;
  1253.     edge_num = j;
  1254.     for ( i = edge_num - 1; i > 0; i-- )
  1255.     { 
  1256.         sorted = 1;
  1257.         for ( j = 0; j < i; j++ )
  1258.             if ( edge[j]->score < edge[j+1]->score )
  1259.             { 
  1260.                 node1 = edge[j];
  1261.                 edge[j] = edge[j+1];
  1262.                 edge[j+1] = node1;
  1263.                 sorted = 0;
  1264.             }
  1265.         if ( sorted )
  1266.             break;
  1267.     }
  1268.     for ( k = 0; k < edge_num; k++ )
  1269.     { 
  1270.         five = edge[k]->host;
  1271.         three = edge[k]->number;
  1272.         orienti = edge[k]->orienti;
  1273.         orientj = edge[k]->orientj;
  1274.         if ( contigs[five].isthree[orienti] &&
  1275.             contigs[three].isfive[orientj] && contigs[five].other != three )
  1276.         { 
  1277.             contigs[five].isthree[orienti] = 0;
  1278.             contigs[three].isfive[orientj] = 0;
  1279.             contigs[five].next[orienti] = three;
  1280.             contigs[five].orient[orienti] = orientj;
  1281.             contigs[five].node[orienti] = edge[k];
  1282.             contigs[three].isthree[(j = 1 - orientj)] = 0;
  1283.             contigs[five].isfive[(i = 1 - orienti)] = 0;
  1284.             contigs[three].next[j] = five;
  1285.             contigs[three].orient[j] = i;
  1286.             contigs[three].node[j] = edge[k];
  1287.             i = contigs[three].other;
  1288.             j = contigs[five].other;
  1289.             contigs[i].other = j;
  1290.             contigs[j].other = i;
  1291.         }
  1292.     }
  1293. }
  1294.  
  1295. RECONCILE(y, pp,nn) overptr y; 
  1296. int  *pp,*nn;
  1297. { 
  1298.     short orienti, orientj;        /* orientation of segments */
  1299.     short orientk, orientd;        /* orientation of segments */
  1300.     int   i, j, k, d, f;            /* row and column indices  */
  1301.     char *ckalloc();            /* space-allocating function */
  1302.     char *A, *B;
  1303.     int  M, N;
  1304.     overptr node1;
  1305.  
  1306.     k = y->host;
  1307.     d = y->number;
  1308.     orientk = y->orienti;
  1309.     orientd = y->orientj;
  1310.     if ( ! contigs[k].isthree[orientk] )
  1311.     { 
  1312.         if ( ! piece[y->ind].kind ) return (0);
  1313.         if ( contigs[d].isthree[orientd] )
  1314.         { 
  1315.             orienti = orientd;
  1316.             i = d;
  1317.             orientj = contigs[k].orient[orientk];
  1318.             j = contigs[k].next[orientk];
  1319.         }
  1320.         else
  1321.             return (0);
  1322.     }
  1323.     else
  1324.         if ( ! contigs[k].isfive[orientk] )
  1325.         { 
  1326.             if ( ! piece[y->ind].kind ) return (0);
  1327.             if ( contigs[d].isfive[orientd] )
  1328.             { 
  1329.                 orienti = contigs[k].orient[1-orientk];
  1330.                 orienti = 1 - orienti;
  1331.                 i = contigs[k].next[1-orientk];
  1332.                 orientj = orientd;
  1333.                 j = d;
  1334.             }
  1335.             else
  1336.                 return (0);
  1337.         }
  1338.         else
  1339.             return (0);
  1340.     A = orienti ? segment[i].seq : segment[i].rev;
  1341.     M = segment[i].length;
  1342.     B = orientj ? segment[j].seq : segment[j].rev;
  1343.     N = segment[j].length;
  1344.     SCORE = 0;
  1345.     big_pass(A,B,M,N,orienti,orientj);
  1346.     if ( SCORE > CUTOFF && ENDI - STARI > over_len && ENDI == M && STARJ == 0 )
  1347.     { 
  1348.         node1 = ( overptr ) ckalloc( (int ) sizeof(over));
  1349.         node1->kind = 1;
  1350.         node1->host = i;
  1351.         node1->number = j;
  1352.         node1->stari = ++STARI;
  1353.         node1->endi = ENDI;
  1354.         node1->orienti = orienti;
  1355.         node1->starj = ++STARJ;
  1356.         node1->endj = ENDJ;
  1357.         node1->orientj = orientj;
  1358.         node1->score = SCORE;
  1359.         piece[*pp].kind = 1;
  1360.         if ( i == d )
  1361.         { 
  1362.             node1->ind = *pp;
  1363.             node1->next = piece[y->ind].list;
  1364.             piece[y->ind].list = node1;
  1365.             piece[*pp].id = j;
  1366.             piece[*pp].list = 0;
  1367.         }
  1368.         else
  1369.         { 
  1370.             node1->ind = y->ind;
  1371.             piece[*pp].list = node1;
  1372.             node1->next = 0;
  1373.             piece[*pp].id = i;
  1374.         }
  1375.         (*nn)++;
  1376.         (*pp)++;
  1377.         f = contigs[k].other;
  1378.         if ( ! contigs[k].isthree[orientk] )
  1379.         { 
  1380.             contigs[j].isfive[orientj] = 1;
  1381.             contigs[j].isthree[1 - orientj] = 1;
  1382.             contigs[k].isthree[orientk] = 1;
  1383.             contigs[k].isfive[1 - orientk] = 1;
  1384.             contigs[f].other = j;
  1385.             contigs[j].other = f;
  1386.         }
  1387.         else
  1388.         { 
  1389.             contigs[i].isthree[orienti] = 1;
  1390.             contigs[i].isfive[1 - orienti] = 1;
  1391.             contigs[k].isfive[orientk] = 1;
  1392.             contigs[k].isthree[1 - orientk] = 1;
  1393.             contigs[f].other = i;
  1394.             contigs[i].other = f;
  1395.         }
  1396.         contigs[k].other = k;
  1397.         return (1);
  1398.     }
  1399.     return (0);
  1400. }
  1401.  
  1402. /* Construct a tree of overlapping-containment segments */
  1403. FORM_TREE()
  1404. { 
  1405.     register  int   i, j, k;    /* index variables */
  1406.     char *ckalloc();        /* space-allocating function */
  1407.     overptr  node1;        /* temporary pointer */
  1408.     short   orient;        /* orientation of segment */
  1409.     int   group;            /* serial number of contigs  */
  1410.     char *A, *B;            /* pointers to segment sequences */
  1411.     int  stari, endi, starj, endj;/* positions where alignment begins */
  1412.     int  M, N;            /* lengths of segment sequences */
  1413.     int   count;            /* temporary variables */
  1414.  
  1415.     mtree = ( struct TTREE * ) ckalloc( seg_num * sizeof(struct TTREE));
  1416.     for ( i = 0; i < seg_num; i++ )
  1417.     { 
  1418.         mtree[i].head = 0;
  1419.         mtree[i].next = mtree[i].child = mtree[i].brother = -1;
  1420.     }
  1421.     for ( group = 0, i = 0; i < seg_num; i++ )
  1422.         if ( segment[i].kind && contigs[i].group < 0 &&
  1423.             ( contigs[i].isfive[1] || contigs[i].isfive[0] ) )
  1424.         { 
  1425.             orient = contigs[i].isfive[1] ? 1 : 0;
  1426.             mtree[i].head = 1;
  1427.             for ( j = i; ;  )
  1428.             { 
  1429.                 contigs[j].group = group;
  1430.                 mtree[j].orient = orient;
  1431.                 SORT(j, orient);
  1432.                 if ( contigs[j].isthree[orient] )
  1433.                     break;
  1434.                 else
  1435.                 { 
  1436.                     k = contigs[j].next[orient];
  1437.                     node1 = contigs[j].node[orient];
  1438.                     if ( j == node1->host )
  1439.                     { 
  1440.                         stari = node1->stari;
  1441.                         endi  = node1->endi;
  1442.                         starj = node1->starj;
  1443.                         endj  = node1->endj;
  1444.                         A = node1->orienti ? segment[j].seq : segment[j].rev;
  1445.                         B = node1->orientj ? segment[k].seq : segment[k].rev;
  1446.                     }
  1447.                     else
  1448.                     { 
  1449.                         M = segment[j].length;
  1450.                         stari = M + 1 - node1->endj;
  1451.                         endi = M + 1 - node1->starj;
  1452.                         N = segment[k].length;
  1453.                         starj = N + 1 - node1->endi;
  1454.                         endj = N + 1 - node1->stari;
  1455.                         A = node1->orientj ? segment[j].rev : segment[j].seq;
  1456.                         B = node1->orienti ? segment[k].rev : segment[k].seq;
  1457.                     }
  1458.                     M = endi - stari + 1;
  1459.                     N = endj - starj + 1;
  1460.                     sapp = S;
  1461.                     last = 0;
  1462.                     al_len = no_mat = no_mis = 0;
  1463.                     (void) diff(&A[stari]-1, &B[starj]-1,M,N,q,q);
  1464.                     count = ( (N = sapp - S) + 1 ) * sizeof(int );
  1465.                     mtree[k].script = ( int  * ) ckalloc( count );
  1466.                     for ( M = 0; M < N; M++)
  1467.                         mtree[k].script[M] = S[M];
  1468.                     mtree[k].size = N;
  1469.                     mtree[k].begin = stari;
  1470.                     mtree[j].next = k;
  1471.                     orient = contigs[j].orient[orient];
  1472.                     j = k;
  1473.                 }
  1474.             }
  1475.             group++;
  1476.         }
  1477. }
  1478.  
  1479. /* Sort the children of each node by the `begin' field */
  1480. SORT(seg, ort) int  seg; 
  1481. short ort;
  1482. { 
  1483.     register  int   i, j, k;    /* index variables */
  1484.     char *ckalloc();        /* space-allocating function */
  1485.     overptr  node1;        /* temporary pointer */
  1486.     short   orient;        /* orientation of segment */
  1487.     char *A, *B;            /* pointers to segment sequences */
  1488.     int  stari, endi, starj, endj;/* positions where alignment begins */
  1489.     int  M, N;            /* lengths of segment sequences */
  1490.     int   count;            /* temporary variables */
  1491.  
  1492.     for ( j = contigs[seg].child; j != -1; j = contigs[j].brother )
  1493.     { 
  1494.         node1 = contigs[j].node[1];
  1495.         if ( ort == node1->orientj )
  1496.         { 
  1497.             stari = node1->starj;
  1498.             endi  = node1->endj;
  1499.             starj = node1->stari;
  1500.             endj  = node1->endi;
  1501.             A = node1->orientj ? segment[seg].seq : segment[seg].rev;
  1502.             B = node1->orienti ? segment[j].seq : segment[j].rev;
  1503.             orient = node1->orienti;
  1504.         }
  1505.         else
  1506.         { 
  1507.             M = segment[seg].length;
  1508.             stari = M + 1 - node1->endj;
  1509.             endi = M + 1 - node1->starj;
  1510.             N = segment[j].length;
  1511.             starj = N + 1 - node1->endi;
  1512.             endj = N + 1 - node1->stari;
  1513.             A = node1->orientj ? segment[seg].rev : segment[seg].seq;
  1514.             B = node1->orienti ? segment[j].rev : segment[j].seq;
  1515.             orient = 1 -  node1->orienti;
  1516.         }
  1517.         M = endi - stari + 1;
  1518.         N = endj - starj + 1;
  1519.         sapp = S;
  1520.         last = 0;
  1521.         al_len = no_mat = no_mis = 0;
  1522.         (void) diff(&A[stari]-1, &B[starj]-1,M,N,q,q);
  1523.         count = ( (M = sapp - S ) + 1 ) * sizeof(int );
  1524.         mtree[j].script = ( int  * ) ckalloc( count );
  1525.         for ( k = 0; k < M; k++)
  1526.             mtree[j].script[k] = S[k];
  1527.         mtree[j].size = M;
  1528.         mtree[j].begin = stari;
  1529.         mtree[j].orient = orient;
  1530.         if ( mtree[seg].child == -1 )
  1531.             mtree[seg].child = j;
  1532.         else
  1533.         { 
  1534.             i = mtree[seg].child;
  1535.             if ( mtree[i].begin >= stari )
  1536.             { 
  1537.                 mtree[j].brother = i;
  1538.                 mtree[seg].child = j;
  1539.             }
  1540.             else
  1541.             { 
  1542.                 M = mtree[i].brother;
  1543.                 for ( ; M != -1; i = M, M = mtree[M].brother )
  1544.                     if ( mtree[M].begin >= stari ) break;
  1545.                 mtree[j].brother = M;
  1546.                 mtree[i].brother = j;
  1547.             }
  1548.         }
  1549.         SORT(j, orient);
  1550.     }
  1551. }
  1552.  
  1553.  
  1554. /* Display the alignments of segments */
  1555. SHOW()
  1556. { 
  1557.     register  int   i, j, k;    /* index variables */
  1558.     char *ckalloc();        /* space-allocating function */
  1559.     int   n;            /* number of working segments */
  1560.     int   limit;            /* number of slots in work */
  1561.     int   col;            /* number of output columns prepared */
  1562.     short done;            /* tells if current group is done */
  1563.     rowptr root;            /* pointer to root of op tree */
  1564.     int   sym[6];            /* occurrance counts for six chars */
  1565.     char  c;            /* temp variable */
  1566.     rowptr t, w, yy;        /* temp pointer */
  1567.     int    x;            /* temp variables */
  1568.     int   group;            /* Contigs number */
  1569.     char  conlit[20], *a;        /* String form of contig number */
  1570.     char  *spt;            /* pointer to the start of consensus */
  1571.  
  1572.     work = ( rowptr * ) ckalloc( seg_num * sizeof(rowptr));
  1573.     group = 0;
  1574.     yy = 0;
  1575.     for ( j = 0; j < 6; j++ )
  1576.         sym[j] = 0;
  1577.     n = limit = col = 0;
  1578.     for ( i = 0; i < seg_num; i++ )
  1579.         if ( mtree[i].head )
  1580.         { 
  1581.             (void) sprintf(conlit, "\n>Contig-%d\n", group);
  1582.             for ( a = conlit; *a; )
  1583.                 *allconpt++ = *a++;
  1584.             /* Mod by S.S.
  1585.      (void) printf("\n#Contig %d\n\n", group++);
  1586. */
  1587.             group++;
  1588.             done = 0;
  1589.             ENTER(&limit, &n, i, col, yy);
  1590.             root = work[0];
  1591.             spt = allconpt;
  1592.             while ( ! done )
  1593.             { 
  1594.                 for ( j = 0; j < n; j++ )    /* get segments into work */
  1595.                 { 
  1596.                     t = work[j];
  1597.                     k = t->id;
  1598.                     if ((x = mtree[k].next) != -1 && mtree[x].begin == t->loc)
  1599.                     { 
  1600.                         ENTER(&limit, &n, x, col, t);
  1601.                         mtree[k].next = -1;
  1602.                     }
  1603.                     for ( x = mtree[k].child; x != -1; x = mtree[x].brother )
  1604.                         if ( mtree[x].begin == t->loc )
  1605.                         { 
  1606.                             ENTER(&limit, &n, x, col, t);
  1607.                             mtree[k].child = mtree[x].brother;
  1608.                         }
  1609.                         else
  1610.                             break;
  1611.                 }
  1612.                 COLUMN(root);        /* determine next column */
  1613.                 root->c = root->kind;
  1614.                 for ( t = head; t != 0; t = t->link )
  1615.                     t->c = t->kind;
  1616.                 for ( j = 0; j < n; j++ )
  1617.                 { 
  1618.                     t = work[j];
  1619.                     if ( t->done )
  1620.                         *t->a++ = ' ';
  1621.                     else
  1622.                     { 
  1623.                         if ( t->c == 'L' )
  1624.                         { 
  1625.                             if ( t->loc == 1 )
  1626.                                 t->offset = allconpt - spt;
  1627.                             c = *t->a++ = t->seq[t->loc++];
  1628.                         }
  1629.                         else
  1630.                             if ( t->loc > 1 )
  1631.                                 c = *t->a++ = '-';
  1632.                             else
  1633.                                 c = *t->a++ = ' ';
  1634.                         if ( c != ' ' )
  1635.                             if ( c == '-' )
  1636.                                 sym[5] += 1;
  1637.                             else
  1638.                                 sym[vert[c]] += 1;
  1639.                         t->c = ' ';
  1640.                     }
  1641.                 }
  1642.                 /* determine consensus char */
  1643.                 k = sym[0] + sym[1] + sym[2] + sym[3] + sym[4];
  1644.                 if ( k < sym[5] )
  1645.                     *allconpt++ = '-';
  1646.                 else
  1647.                     if ( sym[0] == sym[1] && sym[1] == sym[2] &&
  1648.                         sym[2] == sym[3] )
  1649.                         *allconpt++ = 'N';
  1650.                     else
  1651.                     { 
  1652.                         k = sym[0];
  1653.                         c = 'A';
  1654.                         if ( k < sym[1] ) { 
  1655.                             k = sym[1]; 
  1656.                             c = 'C'; 
  1657.                         }
  1658.                         if ( k < sym[2] ) { 
  1659.                             k = sym[2]; 
  1660.                             c = 'G'; 
  1661.                         }
  1662.                         if ( k < sym[3] ) c = 'T';
  1663.                         *allconpt++ = c;
  1664.                     }
  1665.                 for ( j = 0; j < 6; j++ )
  1666.                     sym[j] = 0;
  1667.                 for ( t = head; t != 0; t = t->link )
  1668.                 { 
  1669.                     NEXTOP(t);
  1670.                     if ( t->done )    /* delete it from op tree */
  1671.                     { 
  1672.                         w = t->father;
  1673.                         if ( w->child->id == t->id )
  1674.                             w->child = t->brother;
  1675.                         else
  1676.                         { 
  1677.                             w = w->child;
  1678.                             for ( ; w->brother->id != t->id; w = w->brother )
  1679.                                 ;
  1680.                             w->brother = t->brother;
  1681.                         }
  1682.                     }
  1683.                 }
  1684.                 if ( root->loc > root->length )    /* check root node */
  1685.                 { 
  1686.                     root->done = 1;
  1687.                     if ( (w = root->child) != 0 )
  1688.                     { 
  1689.                         w->father = 0;
  1690.                         root = w;
  1691.                     }
  1692.                     else
  1693.                         done = 1;
  1694.                 }
  1695.                 col++;
  1696.                 if ( col == LINELEN || done )    /* output */
  1697.                 { 
  1698.                     col = 0;
  1699.                     for ( j = 0; j < n; j++ )
  1700.                     { 
  1701.                         t = work[j];
  1702.                         if ( t->done )
  1703.  
  1704. /*
  1705.         Mod by S.S.
  1706.               { (void) printf("#");
  1707.             for ( a = t->name; *a; a++ )
  1708.               (void) printf("%c", *a);
  1709. */
  1710.  
  1711.             {  /* dgg output == fasta compatible */ 
  1712.             short ii, isrev;
  1713.             isrev= (t->name[strlen(t->name)] == '-');
  1714.             
  1715.             (void) fprintf(fout,">cap_");
  1716.             for ( a = t->name; *a; a++ ) (void) fprintf(fout,"%c", *a); 
  1717.             (void) fprintf(fout,"\n");
  1718.             
  1719.             /* need to print offset (.) and reverse as needed */    
  1720.             for ( ii = 0, k=0 ; ii < t->offset; ii++ ) {
  1721.                 k++;
  1722.                 (void) fprintf(fout,"%c", '.'); 
  1723.                 if ( k == LINELEN ) { 
  1724.                     (void)  fprintf(fout,"\n");
  1725.                     k = 0;
  1726.                     }            
  1727.                 }
  1728.         
  1729.             if (isrev) {
  1730.                 for ( a = t->a-1 ; a != t->line-1; a-- )
  1731.                     if ( *a != ' ' ) { 
  1732.                         k++;
  1733.                         (void) fprintf(fout,"%c", nuccomp(*a));
  1734.                         if ( k == LINELEN ) { 
  1735.                             (void)  fprintf(fout,"\n");
  1736.                             k = 0;
  1737.                             }            
  1738.                     }
  1739.                 }
  1740.             else {
  1741.                 for ( /*k = 0,*/ a = t->line ; a != t->a; a++ )
  1742.                     if ( *a != ' ' ) { 
  1743.                         k++;
  1744.                         (void)  fprintf(fout,"%c", *a);
  1745.                         if ( k == LINELEN ) { 
  1746.                             (void)  fprintf(fout,"\n");
  1747.                             k = 0;
  1748.                             }
  1749.                     }
  1750.                 }
  1751.             if (k) (void)  fprintf(fout,"\n");
  1752.             }
  1753.             
  1754. /***** s.s.                
  1755.                         {
  1756.                             int jj;
  1757.                             (void) printf("{\nname    ");
  1758.                             for(jj=0;jj<strlen(t->name)-1;jj++)
  1759.                                 (void) printf("%c", t->name[jj]);
  1760.                             printf("\nstrandedness    %c\n",
  1761.                                 t->name[strlen(t->name)] == '+'? '1':'2');
  1762.  
  1763.                             printf("offset    %d\ntype  DNA\ngroup-ID    %d\nsequence    \"\n",
  1764.                                 t->offset,group);
  1765.                             for ( k = 0, a = t->line ; a != t->a; a++ )
  1766.                                 if ( *a != ' ' )
  1767.                                 { 
  1768.                                     k++;
  1769.                                     (void)  printf("%c", *a);
  1770.                                     if ( k == LINELEN )
  1771.                                     { 
  1772.                                         (void)  printf("\n");
  1773.                                         k = 0;
  1774.                                     }
  1775.                                 }
  1776.  
  1777.                             (void)  printf("\"\n}\n");
  1778.                         }
  1779. ******/
  1780.                         if ( t->linesize - (t->a - t->line) < LINELEN + 3 )
  1781.                             ALOC_SEQ(t);
  1782.                     }
  1783.                     if ( !done )
  1784.                     { 
  1785.                         for ( k = j = n - 1; j >= 0; j-- )
  1786.                             if ( work[j]->done )
  1787.                             { 
  1788.                                 t = work[j];
  1789.                                 for ( x = j; x < k; x++ )
  1790.                                     work[x] = work[x+1];
  1791.                                 work[k--] = t;
  1792.                             }
  1793.                         n = k + 1;
  1794.                     }
  1795.                     else
  1796.                         n = 0;
  1797.                 }
  1798.             }
  1799.         }
  1800.  
  1801.     *allconpt = 0; /* dgg, term string */
  1802.   /* (void)    fprintf(fout,"\n>Consensus\n"); -- dgg added back */
  1803.     for ( a = allconsen, k=0 ; *a; a++ ) {
  1804.         k++;
  1805.         (void) fprintf(fout,"%c", *a); 
  1806.         if ( k == LINELEN ) { 
  1807.             (void) fprintf(fout,"\n");
  1808.             k = 0;
  1809.             }            
  1810.         }
  1811.     fprintf(fout,"\n\n");
  1812. }
  1813.  
  1814. /* allocate more space for output fragment */
  1815. ALOC_SEQ(t) rowptr t;
  1816. { 
  1817.     char  *start, *end, *p;
  1818.     t->linesize *= 2;
  1819.     start = t->line;
  1820.     end = t->a;
  1821.     t->line = ( char * ) ckalloc( t->linesize * sizeof(char));
  1822.     for ( t->a = t->line, p  = start ; p != end; )
  1823.         *t->a++ = *p++;
  1824.     free(start);
  1825. }
  1826.  
  1827. /* enter a segment into working set */
  1828. ENTER(b, d, id, pos, par) int  *b, *d, id, pos; 
  1829. rowptr par;
  1830. { 
  1831.     int  i;
  1832.     char *ckalloc();        /* space-allocating function */
  1833.     rowptr  t;
  1834.  
  1835.     if ( *b <= *d )
  1836.     { 
  1837.         work[*b] = ( rowptr ) ckalloc( (int ) sizeof(row));
  1838.         work[*b]->line = ( char * ) ckalloc( SEQLEN * sizeof(char));
  1839.         work[*b]->linesize = SEQLEN;
  1840.         *b += 1;
  1841.     }
  1842.     t = work[*d];
  1843.     *d += 1;
  1844.     t->a = t->line;
  1845.     for ( i = 0; i < pos; i++ )
  1846.         *t->a++ = ' ';
  1847.     t->c = ' ';
  1848.     t->seq = mtree[id].orient ? segment[id].seq : segment[id].rev;
  1849.     t->length = segment[id].length;
  1850.     t->id = id;
  1851.     if ( par != 0 )
  1852.     { 
  1853.         t->s = mtree[id].script;
  1854.         t->size = mtree[id].size;
  1855.     }
  1856.     t->op = 0;
  1857.     for ( i = 1; i <= segment[id].len && i <= NAMELEN; i++ )
  1858.         t->name[i-1] = segment[id].name[i];
  1859.     if ( mtree[id].orient )
  1860.         t->name[i-1] = '+';
  1861.     else
  1862.         t->name[i-1] = '-';
  1863.     t->name[i] = '\0';
  1864.     t->done = 0;
  1865.     t->loc = 1;
  1866.     t->child = 0;
  1867.     t->father = par;
  1868.     if ( par != 0 )
  1869.     { 
  1870.         t->brother = par->child;
  1871.         par->child = t;
  1872.         NEXTOP(t);
  1873.     }
  1874. }
  1875.  
  1876. /* get the next operation */
  1877. NEXTOP(t) rowptr  t;
  1878. {    
  1879.     if ( t->size || t->op )
  1880.         if ( t->op == 0 && *t->s == 0 )
  1881.         { 
  1882.             t->op = *t->s++;
  1883.             t->size--;
  1884.             t->up = 'L';
  1885.             t->dw = 'L';
  1886.         }
  1887.         else
  1888.         { 
  1889.             if ( t->op == 0 )
  1890.             { 
  1891.                 t->op = *t->s++;
  1892.                 t->size--;
  1893.             }
  1894.             if ( t->op > 0 )
  1895.             { 
  1896.                 t->up = '-';
  1897.                 t->dw = 'L';
  1898.                 t->op--;
  1899.             }
  1900.             else
  1901.             { 
  1902.                 t->up = 'L';
  1903.                 t->dw = '-';
  1904.                 t->op++;
  1905.             }
  1906.         }
  1907.     else
  1908.         if ( t->loc > t->length )
  1909.             t->done = 1;
  1910. }
  1911.  
  1912. COLUMN(x) rowptr x;
  1913. { 
  1914.     rowptr  y;
  1915.     rowptr  start, end;        /* first and last nodes for subtree */
  1916.  
  1917.     if ( x->child == 0 )
  1918.     { 
  1919.         head = tail = 0;
  1920.         x->kind = 'L';
  1921.     }
  1922.     else
  1923.     { 
  1924.         start = end = 0;
  1925.         x->kind = 'L';
  1926.         for ( y = x->child; y != 0; y = y->brother )
  1927.         { 
  1928.             COLUMN(y);
  1929.             if ( x->kind == y->up )
  1930.                 if ( y->kind == y->dw )
  1931.                 { 
  1932.                     if ( head == 0 )
  1933.                     { 
  1934.                         y->link = 0;
  1935.                         head = tail = y;
  1936.                     }
  1937.                     else
  1938.                     { 
  1939.                         y->link = head;
  1940.                         head = y;
  1941.                     }
  1942.                     if ( end == 0 )
  1943.                         start = head;
  1944.                     else
  1945.                         end->link = head;
  1946.                     end = tail;
  1947.                 }
  1948.                 else
  1949.                     if ( y->kind == '-' )
  1950.                     { 
  1951.                         start = head;
  1952.                         end = tail;
  1953.                         x->kind = '-';
  1954.                     }
  1955.                     else
  1956.                     { 
  1957.                         y->link = 0;
  1958.                         y->kind = '-';
  1959.                         if ( end == 0 )
  1960.                             start = end = y;
  1961.                         else
  1962.                         { 
  1963.                             end->link = y;
  1964.                             end = y;
  1965.                         }
  1966.                     }
  1967.             else
  1968.                 if ( y->kind == y->dw )
  1969.                     if ( x->kind == '-' )
  1970.                         ;
  1971.                     else
  1972.                     { 
  1973.                         if ( head == 0 )
  1974.                         { 
  1975.                             y->link = 0;
  1976.                             head = tail = y;
  1977.                         }
  1978.                         else
  1979.                         { 
  1980.                             y->link = head;
  1981.                             head = y;
  1982.                         }
  1983.                         start = head;
  1984.                         end = tail;
  1985.                         x->kind = '-';
  1986.                     }
  1987.                 else
  1988.                     if ( x->kind == '-' )
  1989.                         if ( y->kind == '-' )
  1990.                         { 
  1991.                             if ( end == 0 )
  1992.                             { 
  1993.                                 start = head;
  1994.                                 end = tail;
  1995.                             }
  1996.                             else
  1997.                                 if ( head == 0 )
  1998. /* code folded from here */
  1999.     ;
  2000. /* unfolding */
  2001.                                 else
  2002.                                 { 
  2003. /* code folded from here */
  2004.     end->link = head;
  2005.     end = tail;
  2006. /* unfolding */
  2007.                                 }
  2008.                         }
  2009.                         else
  2010.                             ;
  2011.                     else
  2012.                     { 
  2013.                         start = head;
  2014.                         end = tail;
  2015.                         x->kind = '-';
  2016.                     }
  2017.         }
  2018.         head = start;
  2019.         tail = end;
  2020.     }
  2021. }
  2022.  
  2023. /* Display a summary of contigs */
  2024. GRAPH()
  2025. { 
  2026.     int   i, j, k;        /* index variables */
  2027.     int   group;            /* serial number of contigs  */
  2028.     char name[NAMELEN+2];        /* name of segment */
  2029.     char  *t;            /* temp var */
  2030.     int   length;            /* length of name */
  2031.  
  2032.     (void)    printf("\nOVERLAPS            CONTAINMENTS\n\n");
  2033.     group = 1;
  2034.     for ( i = 0; i < seg_num; i++ )
  2035.         if ( mtree[i].head )
  2036.         { 
  2037.             (void) printf("******************* Contig %d ********************\n",
  2038.                 group++ );
  2039.             for ( j = i; j != -1; j = mtree[j].next )
  2040.             { 
  2041.                 length = segment[j].len;
  2042.                 t = segment[j].name + 1;
  2043.                 for ( k = 0; k < length && k < NAMELEN; k++ )
  2044.                     name[k] = *t++;
  2045.                 if ( mtree[j].orient )
  2046.                     name[k] = '+';
  2047.                 else
  2048.                     name[k] = '-';
  2049.                 name[k+1] = '\0';
  2050.                 (void) printf("%s\n", name);
  2051.                 CONTAIN(mtree[j].child, name);
  2052.             }
  2053.         }
  2054. }
  2055.  
  2056. CONTAIN(id, f) int  id; 
  2057. char *f;
  2058. { 
  2059.     int   k;            /* index variable */
  2060.     char name[NAMELEN+2];        /* name of segment */
  2061.     char  *t;            /* temp var */
  2062.     int   length;            /* length of name */
  2063.  
  2064.     if ( id != -1 )
  2065.     { 
  2066.         length = segment[id].len;
  2067.         t = segment[id].name + 1;
  2068.         for ( k = 0; k < length && k < NAMELEN; k++ )
  2069.             name[k] = *t++;
  2070.         if ( mtree[id].orient )
  2071.             name[k] = '+';
  2072.         else
  2073.             name[k] = '-';
  2074.         name[k+1] = '\0';
  2075.         (void) printf("                    %s is in %s\n", name,f);
  2076.         CONTAIN(mtree[id].child, name);
  2077.         CONTAIN(mtree[id].brother, f);
  2078.     }
  2079. }
  2080.  
  2081. big_pass(A,B,M,N,orienti,orientj) char A[],B[]; 
  2082. int  M,N; 
  2083. short orienti, orientj;
  2084. { 
  2085.     register  int   i, j;            /* row and column indices */
  2086.     register  int   c;            /* best score at current point */
  2087.     register  int   f;            /* best score ending with insertion */
  2088.     register  int   d;            /* best score ending with deletion */
  2089.     register  int   p;            /* best score at (i-1, j-1) */
  2090.     register  int   ci;            /* end-point associated with c */
  2091.  
  2092.     register  int   di;            /* end-point associated with d */
  2093.     register  int   fi;            /* end-point associated with f */
  2094.     register  int   pi;            /* end-point associated with p */
  2095.     char   *va;                /* pointer to varray(A[i], B[j]) */
  2096.     int   x1, x2;        /* regions of A before x1 or after x2 are lightly penalized */
  2097.     int   y1, y2;        /* regions of B before y1 or after y2 are lightly penalized */
  2098.     short  heavy;        /* 1 = heavy penalty */
  2099.     int   ex, gx;        /* current gap penalty scores */
  2100.  
  2101.     /* determine x1, x2, y1, y2 */
  2102.     if ( POS5 >= POS3 )
  2103.         fatal("The value for POS5 must be less than the value for POS3");
  2104.     if ( orienti )
  2105.     { 
  2106.         x1 = POS5 >= M ? 1 : POS5;
  2107.         x2 = POS3 >= M ? M : POS3;
  2108.     }
  2109.     else
  2110.     { 
  2111.         x1 = POS3 >= M ? 1 : M - POS3 + 1;
  2112.         x2 = POS5 >= M ? M : M - POS5 + 1;
  2113.     }
  2114.     if ( orientj )
  2115.     { 
  2116.         y1 = POS5 >= N ? 1 : POS5;
  2117.         y2 = POS3 >= N ? N : POS3;
  2118.     }
  2119.     else
  2120.     { 
  2121.         y1 = POS3 >= N ? 1 : N - POS3 + 1;
  2122.         y2 = POS5 >= N ? N : N - POS5 + 1;
  2123.     }
  2124.     if ( x1 + 1 <= x2 ) x1++;
  2125.     if ( y1 + 1 <= y2 ) y1++;
  2126.     heavy = 0;
  2127.  
  2128.     /* Compute the matrix.
  2129.        CC : the scores of the current row
  2130.        RR : the starting point that leads to score CC
  2131.        DD : the scores of the current row, ending with deletion
  2132.        SS : the starting point that leads to score DD        */
  2133.     /* Initialize the 0 th row */
  2134.     for ( j = 1; j <= N ; j++ )
  2135.     {  
  2136.         CC[j] = 0;
  2137.         DD[j] = - (q);
  2138.         RR[j] = SS[j] = -j;
  2139.     }
  2140.     for ( i = 1; i <= M; i++)
  2141.     {  
  2142.         if ( i == x1 ) heavy = 1 - heavy;
  2143.         if ( i == x2 ) heavy = 1 - heavy;
  2144.         ex = r1;
  2145.         gx = qr1;
  2146.         va = v1[A[i]];
  2147.         c = 0;                /* Initialize column 0 */
  2148.         f = - (q);
  2149.         ci = fi = i;
  2150.         p = 0;
  2151.         pi = i - 1;
  2152.         for ( j = 1 ; j <= N ; j++ )
  2153.         {  
  2154.             if ( j == y1 )
  2155.             { 
  2156.                 if ( heavy )
  2157.                 { 
  2158.                     ex = r;
  2159.                     gx = qr;
  2160.                     /*
  2161. S.S.
  2162.                     va = varray[A[i]];
  2163. */
  2164.                     va = varray[A[i]];
  2165.                 }
  2166.             }
  2167.             if ( j == y2 )
  2168.             { 
  2169.                 if ( heavy )
  2170.                 { 
  2171.                     ex = r1;
  2172.                     gx = qr1;
  2173.                     va = v1[A[i]];
  2174.                 }
  2175.             }
  2176.             if ( ( f = f - ex ) < ( c = c - gx ) )
  2177.             { 
  2178.                 f = c; 
  2179.                 fi = ci; 
  2180.             }
  2181.             di = SS[j];
  2182.             if ( ( d = DD[j] - ex ) < ( c = CC[j] - gx ) )
  2183.             { 
  2184.                 d = c; 
  2185.                 di = RR[j]; 
  2186.             }
  2187.             c = p+va[B[j]];        /* diagonal */
  2188.             ci = pi;
  2189.             if ( c < d )
  2190.             { 
  2191.                 c = d; 
  2192.                 ci = di; 
  2193.             }
  2194.             if ( c < f )
  2195.             { 
  2196.                 c = f; 
  2197.                 ci = fi; 
  2198.             }
  2199.             p = CC[j];
  2200.             CC[j] = c;
  2201.             pi = RR[j];
  2202.             RR[j] = ci;
  2203.             DD[j] = d;
  2204.             SS[j] = di;
  2205.             if ( ( j == N || i == M ) &&  c > SCORE )
  2206.             { 
  2207.                 SCORE = c;
  2208.                 ENDI = i;
  2209.                 ENDJ = j;
  2210.                 STARI = ci;
  2211.             }
  2212.         }
  2213.     }
  2214.     if ( SCORE )
  2215.         if ( STARI < 0 )
  2216.         { 
  2217.             STARJ = - STARI;
  2218.             STARI = 0;
  2219.         }
  2220.         else
  2221.             STARJ = 0;
  2222. }
  2223.  
  2224. /* diff(A,B,M,N,tb,te) returns the score of an optimum conversion between
  2225.    A[1..M] and B[1..N] that begins(ends) with a delete if tb(te) is zero
  2226.    and appends such a conversion to the current script.                   */
  2227.  
  2228. int  diff(A,B,M,N,tb,te) char *A, *B; 
  2229. int  M, N; 
  2230. int  tb, te;
  2231.  
  2232. { 
  2233.     int    midi, midj, type;    /* Midpoint, type, and cost */
  2234.     int  midc;
  2235.  
  2236.     { 
  2237.         register int    i, j;
  2238.         register int  c, e, d, s;
  2239.         int  t;
  2240.         char    *va;
  2241.         char  *ckalloc();
  2242.  
  2243.         /* Boundary cases: M <= 1 or N == 0 */
  2244.  
  2245.         if (N <= 0)
  2246.         { 
  2247.             if (M > 0) DEL(M)
  2248.                 return - gap(M);
  2249.         }
  2250.         if (M <= 1)
  2251.         { 
  2252.             if (M <= 0)
  2253.             { 
  2254.                 INS(N);
  2255.                 return - gap(N);
  2256.             }
  2257.             if (tb > te) tb = te;
  2258.             midc = - (tb + r + gap(N) );
  2259.             midj = 0;
  2260.             va = varray[A[1]];
  2261.             for (j = 1; j <= N; j++)
  2262.             { 
  2263.                 c = va[B[j]] - ( gap(j-1) + gap(N-j) );
  2264.                 if (c > midc)
  2265.                 { 
  2266.                     midc = c;
  2267.                     midj = j;
  2268.                 }
  2269.             }
  2270.             if (midj == 0)
  2271.             { 
  2272.                 INS(N) DEL(1)             }
  2273.             else
  2274.             { 
  2275.                 if (midj > 1) INS(midj-1)
  2276.                     REP
  2277.                         if ( (A[1]|32) == (B[midj]|32) )
  2278.                         no_mat += 1;
  2279.                     else
  2280.                         no_mis += 1;
  2281.                 if (midj < N) INS(N-midj)
  2282.             }
  2283.             return midc;
  2284.         }
  2285.  
  2286.         /* Divide: Find optimum midpoint (midi,midj) of cost midc */
  2287.  
  2288.         midi = M/2;            /* Forward phase:                          */
  2289.         CC[0] = 0;            /*   Compute C(M/2,k) & D(M/2,k) for all k */
  2290.         t = -q;
  2291.         for (j = 1; j <= N; j++)
  2292.         { 
  2293.             CC[j] = t = t-r;
  2294.             DD[j] = t-q;
  2295.         }
  2296.         t = -tb;
  2297.         for (i = 1; i <= midi; i++)
  2298.         { 
  2299.             s = CC[0];
  2300.             CC[0] = c = t = t-r;
  2301.             e = t-q;
  2302.             va = varray[A[i]];
  2303.             for (j = 1; j <= N; j++)
  2304.             { 
  2305.                 if ((c = c - qr) > (e = e - r)) e = c;
  2306.                 if ((c = CC[j] - qr) > (d = DD[j] - r)) d = c;
  2307.                 c = s+va[B[j]];
  2308.                 if (c < d) c = d;
  2309.                 if (c < e) c = e;
  2310.                 s = CC[j];
  2311.                 CC[j] = c;
  2312.                 DD[j] = d;
  2313.             }
  2314.         }
  2315.         DD[0] = CC[0];
  2316.  
  2317.         RR[N] = 0;            /* Reverse phase:                          */
  2318.         t = -q;            /*   Compute R(M/2,k) & S(M/2,k) for all k */
  2319.         for (j = N-1; j >= 0; j--)
  2320.         { 
  2321.             RR[j] = t = t-r;
  2322.             SS[j] = t-q;
  2323.         }
  2324.         t = -te;
  2325.         for (i = M-1; i >= midi; i--)
  2326.         { 
  2327.             s = RR[N];
  2328.             RR[N] = c = t = t-r;
  2329.             e = t-q;
  2330.             va = varray[A[i+1]];
  2331.             for (j = N-1; j >= 0; j--)
  2332.             { 
  2333.                 if ((c = c - qr) > (e = e - r)) e = c;
  2334.                 if ((c = RR[j] - qr) > (d = SS[j] - r)) d = c;
  2335.                 c =  s+va[B[j+1]];
  2336.                 if (c < d) c = d;
  2337.                 if (c < e) c = e;
  2338.                 s = RR[j];
  2339.                 RR[j] = c;
  2340.                 SS[j] = d;
  2341.             }
  2342.         }
  2343.         SS[N] = RR[N];
  2344.  
  2345.         midc = CC[0]+RR[0];        /* Find optimal midpoint */
  2346.         midj = 0;
  2347.         type = 1;
  2348.         for (j = 0; j <= N; j++)
  2349.             if ((c = CC[j] + RR[j]) >= midc)
  2350.                 if (c > midc || CC[j] != DD[j] && RR[j] == SS[j])
  2351.                 { 
  2352.                     midc = c;
  2353.                     midj = j;
  2354.                 }
  2355.         for (j = N; j >= 0; j--)
  2356.             if ((c = DD[j] + SS[j] + q) > midc)
  2357.             { 
  2358.                 midc = c;
  2359.                 midj = j;
  2360.                 type = 2;
  2361.             }
  2362.     }
  2363.  
  2364.     /* Conquer: recursively around midpoint */
  2365.  
  2366.     if (type == 1)
  2367.     {
  2368.         (void) diff(A,B,midi,midj,tb,q);
  2369.         (void) diff(A+midi,B+midj,M-midi,N-midj,q,te);
  2370.     }
  2371.     else
  2372.     {
  2373.         (void) diff(A,B,midi-1,midj,tb,zero);
  2374.         DEL(2);
  2375.         (void) diff(A+midi+1,B+midj,M-midi-1,N-midj,zero,te);
  2376.     }
  2377.     return midc;
  2378. }
  2379.  
  2380. /* lib.c - library of C procedures. */
  2381.  
  2382. /* fatal - print message and die */
  2383. fatal(msg)
  2384. char *msg;
  2385. {
  2386.     (void) fprintf(stderr, "%s\n", msg);
  2387.     exit(1);
  2388. }
  2389.  
  2390. /* fatalf - format message, print it, and die */
  2391. fatalf(msg, val)
  2392. char *msg, *val;
  2393. {
  2394.     (void) fprintf(stderr, msg, val);
  2395.     (void)    putc('\n', stderr);
  2396.     exit(1);
  2397. }
  2398.  
  2399. /* ckopen - open file; check for success */
  2400. FILE *ckopen(name, mode)
  2401. char *name, *mode;
  2402. {
  2403.     FILE *fopen(), *fp;
  2404.  
  2405.     if ((fp = fopen(name, mode)) == NULL)
  2406.         fatalf("Cannot open %s.", name);
  2407.     return(fp);
  2408. }
  2409.  
  2410. /* ckalloc - allocate space; check for success */
  2411. char *ckalloc(amount)
  2412. int  amount;
  2413. {
  2414.     char *p;
  2415.  
  2416.     if ((p = (char*) malloc( (unsigned) amount)) == NULL)
  2417.         fatal("Ran out of memory.");
  2418.     return(p);
  2419. }
  2420.  
  2421.